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The lack of regularity of geometric variables at the origin is often a source of serious problem 
for spherically symmetric evolution codes in numerical relativity. One usually deals with this by 
restricting the gauge and solving the hamiltonian constraint for the metric. Here we present a 
generic algorithm for dealing with the regularization of the origin that can be used directly on the 
evolution equations and that allows very general gauge choices. Our approach is similar in spirit 
to the one introduced by Arbona and Bona for the particular case of the Bona-Masso formulation. 
However, our algorithm is more general and can be used with a wide variety of evolution systems. 
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I. INTRODUCTION 

When developing spherically symmetric codes in nu- 
merical relativity, the coordinate singularity at the origin 
can be a source of serious problems caused by the lack 
of regularity of the geometric variables there. The prob- 
lem arises because of the presence of terms in the evolu- 
tion equations that go as 1 [r near the origin. Regularity 
of the metric (essentially local flatness) guarantees the 
exact cancellations of such terms at the origin, thus en- 
suring well-behaved solutions. This exact cancellation, 
however, though certainly true for analytical solutions, 
usually fails to hold for numerical solutions. One then 
finds that the 1 /r terms do not cancel and the numerical 
solution becomes ill-behaved near r — 0: it not only fails 
to converge there, but can easily turn out to be violently 
unstable in just a few timesteps. 

The usual way to deal with this problem is to use the 
so-called areal (or radial) gauge, where the radial coor- 
dinate r is chosen in such a way that the proper area of 
spheres of constant r is always Aitr 2 . If, moreover, one 
also choses a vanishing shift vector one ends up in the 
standard polar /areal gauge 0,0 j for which the lapse is 
forced to satisfy a certain ordinary differential equation in 
r. The name polar comes from the fact that for this gauge 
choice there is only one non-zero component of the extrin- 
sic curvature tensor, namely K rr In the polar/areal 
gauge the problem of achieving the exact cancellation of 
the 1/r terms is reduced to imposing the boundary con- 
dition g rr = 1 at r = 0, which can be easily done if one 
solves for g rr from the hamiltonian constraint (which in 
this case is an ordinary differential equation in r) and ig- 
nores its evolution equation. If one does this in vacuum, 
one ends up inevitably with Minkowski spacetime in the 
usual coordinates (one can also recover Schwarzschild by 
working in isotropic coordinates and factoring out the 
conformal factor analytically). Of course, in the presence 
of matter, one can still have truly non-trivial dynamics. 
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The main drawback of the standard approach is that 
the gauge choice has been completely exhausted. In par- 
ticular, the polar/areal gauge can not penetrate apparent 
horizons, since inside an apparent horizon it is impossible 
to keep the areas of spheres fixed without a non-trivial 
shift vector. The polar/areal gauge has, nevertheless, 
been used successfully even in the study of critical col- 
lapse to a black hole, where the presence of the black hole 
is identified by the familiar "collapse of the lapse" even 
if no apparent horizon can be found 0]. 

Still, one would like to have a way of dealing with the 
regularity issue that allows more generic gauge choices to 
be made, either because one is interested in studying the 
region inside an apparent horizon, or because one wants 
to test interesting gauge conditions in the simple case of 
spherical symmetry. Because of this we have developed a 
general regularization technique that can be used directly 
on the Einstein evolution equations. 

Our regularization method is similar in spirit, if not in 
detail, to the one presented by Arbona and Bona in 
The main difference being that the approach of Arbona 
and Bona was tied to the use of the Bona-Masso evolution 
system [U El 0, E 13 j while our algorithm is much more 
general. 

A final point deserves notice. Spherically symmetric 
evolution codes that involve eternal black holes usually 
ignore the regularity problem. For example, one can ex- 
cise the black hole interior and eliminate r — from the 
numerical grid. But even if one does not use excision, 
for a black hole r — is not a regular point but rather a 
compactification of the asymptotic infinity on the other 
side of the Einstein-Rosen bridge. This compactifica- 
tion introduces geometric factors that compensate the 
1/r terms making the equations regular even at r = 0. 
This means that, contrary to what one would naively ex- 
pect, in spherical symmetry it is easier to evolve eternal 
black holes (at least for some time) than it is to evolve 
regular spacetimes. 

This paper is organized as follows. In Sec. [H] we dis- 
cuss the regularity conditions that the metric functions 
must satisfy at the origin of spherical coordinates, and 
we show which terms need to be regularized in the Ein- 
stein equations. Section llTTl describes our regularization 
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algorithm in a generic way. In Sec. HVI we present exam- 
ple of how to regularize some specific formulations of the 
Einstein equations, and show some numerical examples. 
We conclude in section Ivl 



II. REGULARITY CONDITIONS 

We start by writing the general form of the spatial 
metric in spherical symmetry as 



dl 2 



A(r,t)dr 2 +r 2 B(r,t)dQ 2 



(2.1) 



with A and B positive metric functions and dVl 2 the solid 



angle element: dfl — dO + sin 



Notice that we 



have already factored out the r 2 dependency of the an- 
gular metric functions. This has the advantage of making 
explicit the dependency on r of geometric quantities and 
makes the regularization procedure easier. 

As we will deal with the Einstein equations in first 
order form, we will introduce the auxiliary quantities: 



D A := d r In A 



D B := d r \i\B 



(2.2) 



We will also work with the mixed components of the ex- 
trinsic curvature: K A ■= K r r , Kb ■= Kg — Kt. 

There are in fact two different types of regularity con- 
ditions that the variables {A, B, Da, Db,K a , Kb} must 
satisfy at r = 0. The first type of conditions are simple 
those imposed by the requirement that the different vari- 
ables should be well defined at the origin, and imply the 
following behavior for small r: 



A r 


- A° + 


0(r 2 ) , 


(2.3) 


B r- 


- B° + 


0(r 2 ) , 


(2.4) 


D A - 


- 0(r) 


) 


(2.5) 


D B - 


- C(r) 


) 


(2.6) 


K A - 


- 


-0(r 2 ), 


(2.7) 


K B - 


- K%1 


-0(r 2 ), 


(2.8) 



with {A , B°, K A , K B } perhaps functions of time, but 
not of r. These regularity conditions are in fact quite 
easy to implement numerically. For example, one can 
use a finite differencing grid that staggers the origin, and 
then obtain data on the fictitious point at r — — Ar/2 
by demanding for {A, B, Ka, Kb} to be even functions 
at r — and for {Da, Db} to be odd. 

It is the second type of regularity conditions that is 
more troublesome. In order to see the problem, we 
will first write the Arnowitt-Deser-Misner (ADM) equa- 
tions [ToL ITU in first order form for the case of spherical 
symmetry. The evolution equations are (in the case of 
zero shift) 



d t A = -IolAKa , 
d t B = -2aBK B , 
d t D A = -2a[K A D a 



d r K A ] 



(2.9) 
(2.10) 
(2.11) 



d t D B 



8 t K A = - A 
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2a[K B D a + d r K B ] , 

d r (D a +D B )+D 2 a - 
D A D B 



(2.12) 



D a D A 



- AK A (K A + 2Kb) 



(D A -2D B ) 



(2.13) 



d t K E 



a 

2A 



d r D B + D a D B + D 2 B 



D A D L 



I, s 2(A-B) 

- -{D A - 2D a - AD B ) - V 2p ; 
r r B 

+ aK B {K A + 2K B ) 



(2.14) 



where a is the lapse function and D a :— d r \na. The 
hamiltonian and momentum constraints take the form 



d r D B = 



1 

r^B 



(A — B)+ AK B {2K A + Kb) 



1 -(D a -3Db) + ^-^,(2.1 5 ) 



d r K B = [K A -K L 



£b 
2 



(2.16) 



Since {D a , D A , D B } go as r near the origin, terms of 
the type D{ aA B}/r are regular and represent no prob- 
lem. However, we see that both in the hamiltonian con- 
straint and in the evolution equation for Kb there is a 
term of the form (A—B)/r 2 , while in the momentum con- 
straint there is a term of the form (K A — Kb) It. Given 
the behavior of these variables near r = these terms 
would seem to blow up at the origin. The reason why 
this does not in fact happen is that, near the origin, we 
must also ask for the extra regularity conditions 



A-B ~<D{r 2 ) 



that is 



A° = B° 



K A -K B ~ 0{r 2 ) , (2.17) 



K° A = K° B . (2.18) 



It is not difficult to understand where these conditions 
come from. They are just a consequence of the fact that 
space must remain locally flat at r = 0. This local flat- 
ness condition implies that, near r — 0, it must be possi- 
ble to write the metric as 



dl 2 



dR 2 + R 2 dn 2 



(2.19) 



with R a radial coordinate that measures proper distance 
from the origin. A local transformation of coordinates 
from R to r then takes the metric into the form 



dl' 



r~0 



dR 
dr 



(dr 2 + r 2 dn 2 ) 



(2.20) 



r=0 



which implies that A = B° and, since this must hold for 

1 A 



all time, also that K A = K B . 
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It turns out that it is not trivial to implement numeri- 
cally both the symmetry regularity conditions l|2.3|) - (|2.8|) 
and the local flatness regularity conditions (|2.17(1 at the 
same time. The reason for this is that at r — we now 
have three boundary conditions for just two variables: 
both the derivatives of A and B must vanish, plus A and 
B must be equal to each other (and the same thing must 
happen for K A and Kb)- The boundary conditions for 
the exact equations are, of course, also over-determined, 
but in that case the consistency of the equations implies 
that if they are satisfied initially they remain satisfied 
for all time. In the numerical case, however, this is not 
true due to truncation errors, and very rapidly (typically 
within one or two time steps) one of the three bound- 
ary conditions fails to hold. It is easy to convince one- 
self that simply ignoring one condition and imposing the 
other two does not work. If we impose the zero deriva- 
tive condition and ignore the A — B condition, then the 
(A — B) jr 2 term in the evolution equations rapidly be- 
comes singular. On the other hand, if we impose the 
A = B condition and ignore one of the symmetry condi- 
tions, then we introduce an inconsistency with the finite 
difference version of the evolution equations, since for fi- 
nite Ar it is very difficult to guarantee that the difference 
between dtK A and dtKs approaches zero at the origin. 
This inconsistency then very rapidly causes large (and 
non-convergent) gradients to develop near the origin. In 
the following section we will introduce an algorithm that 
successfully regularizes the numerical evolution equations 
near r = 0. 

As a final comment, from the above equations we can 
also very easily see why the polar/areal gauge has no se- 
rious regularity problem. In that gauge we have B = 1 
by construction. If we now impose the boundary condi- 
tion A(r = 0) = 1, and solve for A(r) by integrating the 
hamiltonian constraint with B = 1 and Db = (ignor- 
ing the evolution equations), then the (A — B)/r 2 term 
causes no trouble. 



III. REGULARIZATION ALGORITHM 

In Ref. 0, Arbona and Bona developed a regulariza- 
tion technique for the spherically symmetric version of 
the Bona-Masso (BM) evolution system 0, H H @. 
Their technique is based on redefining the auxiliary dy- 
namical variable V r that is part of the standard BM 
formulation in a way that allows them to absorb the 
(A — B)/r 2 terms and reduces the regularization prob- 
lem to applying the correct boundary condition to V r at 
r = 0. 

Since we are interested in developing a generic regular- 
ization technique, we will start from the ADM equations 
from the previous section. However, we will take the idea 
from the regularization technique of Arbona and Bona of 
introducing an auxiliary variable that will allow us to ab- 
sorb the problematic terms. Adding an auxiliary variable 
seems to us to be the more straightforward way of solving 



the problem of having over-determined boundary condi- 
tions: the extra boundary condition will be imposed on 
the auxiliary variable. We will then define the variable, 



A 



r B 



(3.1) 



Notice that, if the local flatness regularity condi- 
tions (|2~T7|) are satisfied, then the variable A has the 
following behavior at the origin 



A ~ 0(r) , 



(3.2) 



which, as mentioned above, can easily be imposed numer- 
ically using a grid that staggers the origin, and asking for 
A to be odd across r = 0. 

The hamiltonian constraint now becomes 



d r D B = -+AKb{2K a + K b ) 
r 

+ 1 -(D A -3D B ) + ^ 
r 2 



3D% 



(3.3) 



and the evolution equation for Kb becomes 



d t K B 



a 

2A 



d r D B + D a D B + Dj 



D A D l 



-{D A -2D a ~4D B -2\) 
r 

aK B (K A + 2K B ) • 



(3.4) 



As mentioned, the problem terms have now been trans- 
formed into A/r, which will be well behaved as long as 
A is odd at r = 0. The momentum constraint still has 
a term (K A — K B )/r, but this should cause no trouble 
since it does not feed back into the evolution equations 
(one can always multiply the momentum constraint with 
r before evaluating it). Of course, multiples of the mo- 
mentum constraint are typically added to the evolution 
equations in order to build hyperbolic formulations, and 
we will discuss how to deal with that term in such a case 
below. 

There is one other ingredient that needs to be added: 
an evolution equation for A. This can be obtained di- 
rectly from its definition: 



2a A ( K A — Kr 



B 



(3.5) 



The last evolution equation clearly has the dangerous 
(K A — K B )/r term, but this term can be removed with 
the help of the momentum constraint (|2.16|) to find 



2aA 
~B~ 



d r K B -^-{K A -K B ) 



(3.6) 



which is now regular at the origin. 

The regularized first order ADM evolution equations 
are then (|2~§j) - ((2~T3j) . with (|2~T4l) replaced by plus 
the evolution equation for A given by Q3.6p . 
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Having regularized the standard ADM equations, the 
question arises of how to regularize alternative formula- 
tions where multiples of the constraints can be added to 
the evolution equations in a number of ways. Adding 
multiples of the hamiltonian constraint represents no 
problem, as the introduction of the variable A already 
regularized this constraint, as seen in i|3.3[l . The mo- 
mentum constraint, however, is not regularized as it still 
includes the term (Ka — Kb) jr. One could try to play 
the same game as before and introduce yet another vari- 
able to absorb this term. However, we will now show that 
this is not really necessary. 

Let us then consider some arbitrary first order formula- 
tion of the Einstein evolution equations in spherical sym- 
metry that has the generic form 

d f Ui = qi(u,v) , (3.7) 
d t Vi = M({u) d r Vj +pi(u,v) , (3.8) 

where u = (A,B,X) and v = (Da, D b , Ka, K b )- The 
source terms q and p are assumed not to depend on 
derivatives of any of the fields. The formulation might be 
hyperbolic or not, depending on the characteristic struc- 
ture of the matrix M. We will assume that one has ar- 
rived at such a formulation by adding multiples of the 
hamiltonian and momentum constraints to the evolution 
equations for the v's. This means that one can expect 
that the source terms pi will in general contain terms 
proportional to (K a — Kb) jr. We will then rewrite the 
evolution equations for the v, L as 

dtVi = M((u) d r v 3 +p' i (u,v) + ^-(K A - K B ) . (3.9) 

Here we are assuming that the coefficient fi(u) of the 
(K a — K b) i 'r terms depends of the it's, but not on the v's, 
which will typically be the case. Using now equation (|3.5|) 
we find 

d t v t = M((u) d rVj +p' i (u,v) + d t X , (3.10) 



2aA 



which implies 



2aA 



If we now define 



h{u)B 



(3.12) 



The final step is to substitute the spatial derivative of Vj 
for that of v[- to find 

d t v[ = M\{u) drv'j + p'^u, v) - A Fl{u, v) 
+ d ( h[u)B A 
= M\ (u) drv'j + p[{u, v) + A (F[(u, v) - Ff(u, v)) 

(3.14) 

with F[(u,v) = d r (fi(u)B /2a A). Using now the fact 
that 



we finally find 



1 



X+^(D A ~D B ) 



(3.15) 



d t v[ = M 1 i (u)d r v' J +p' l (u,v) + X(F[(u,v)-F l t (u,v)) 
fi( u ) B I,. , A m (3 . 16) 



2aAr 



This last equation is now regular, and has precisely 
the same characteristic structure as the original system. 
What we have done is transform the original evolution 
equations for the variables into evolution equations for 
the new v[ variables for which the principal part terms are 
the same and the source terms are regular. Notice that 
typically only some of the fi(u) will be different from 
zero, so one does not need to transform all variables. 

In the following section we will consider examples of 
how to regularize some specific systems of evolution equa- 
tions. 



IV. EXAMPLES 

The regularized first order ADM equations where de- 
rived in the last section. Here we will examples of numer- 
ical evolutions using two different regularized systems. 

In the numerical simulations we will take as initial data 
Minkowski spacetime in the usual coordinates, so that: 



.4 
D A 
K a 



B = 1 , 
D B = 0, 
K B = 0. 



(4.1) 
(4.2) 
(4.3) 



In order to have a non-trivial evolution, we will chose an 
initial lapse profile of the form: 



a(t = 0) = 1 + r 2 Ce 



(4.4) 



we can transform the last equation into 

d t v[ = Mf{u) d r v 3 +p' l (u,v) - X Ff(u,v) , (3.13) 

with Ff{u,v) = d t (f l (u)B/2aA). Notice that Ff(u,v) 
so defined will involve no spatial derivatives of it's or v's. 



that is, we add a small gaussian contribution to the initial 
Minkowski lapse. We will then evolve the lapse using a 
Bona-Masso (BM) slicing condition so the evolution 
equation for the lapse will be: 



d t a = -a 2 f(a)(KA + 2K B ) . 



(4.5) 



FIG. 1: System I not regularized. The plots show the evo- 
lution of the metric function A at different times. Note the 
spike at r = for t = 1. 



FIG. 2: System I not regularized. The plots show the evolu- 
tion of the lapse function a at different times. Note the spike 
at r = for t = 1. 



In the simulations shown below we have taken the 
gaussian parameters to be: C — 0.001, r — 5.0, a = 1.0. 
We have also restricted ourselves to harmonic slicing, 
that is, f(a) = 1. 



A. System I 

As a first example we will build a hyperbolic system 
starting from the ADM equations and the BM slicing 
condition. We construct this system by using the Hamil- 
tonian and Momentum constraints to remove the terms 
proportional to d r Ds and d r Ks from the evolution equa- 
tions of K A and D a respectively. The resulting system 
can be easily shown to be strongly hyperbolic. 

Figures ^and 121 show the evolution of the radial metric 
A and lapse function a using the system desribed above, 
with no regularization. Note that both plots show a spike 
at r = for times t » 1. 

Next we look at the regularized case. As described 
above, we first introduce the auxilliary variable 



A = 



1 



(4.6) 



Also, since we have used the momentum constraint to 
modify the evolution equation for D a , we will need to 
replace this variable with 



U a := D c 



BX 
~A 



(4.7) 



The set of variables to be evolved is then: 



{a,A,B,U a ,D A ,D B ,K A ,K B ,\} . (4.8) 
The evolution equations for these variables take the form: 
d t a = -a 2 f(K A + 2K B ) , (4.9) 



d t A = -2aAK A , 
d t B = -2aBK B , 

d t U a = -af 8 r K A + a(K A + 2K B ) 



(4.10) 
(4.11) 



A 



U a (/ + af) 
2BX 

D B - 



af(K B - K A ) 



A 



d t D A = -2a 
d t D B = -2a 



(4.12) 
(4.13) 



(4.14) 



o t K A = -—d r U a - —i lU a —J 



^ (/ + «/') 
A (K% K\) 



U 

2 

fBX , , D% 



1 



X + D B 



d t K B 



f(^+D A - Dl 



^^A 



(4.15) 



D A D 



B 



2 

2fBX 



2AKb(K a + 2K b ) 



fBX\ 
1 

1 



.4 



- AD B - 2X 



2aA „ aAD B , 

d t X = — 8 r K B + -^{K B - K A ) 



B 



B 



D A - 2U a 
(4.16) 
(4.17) 



with /' := df I da. 

Figures |2| and 0] shown again the evolution of A and a 
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FIG. 3: System I regularized. The plots show the evolution 
of the metric function A at different times. 




FIG. 4: System I regularized. The plots show the evolution 
of the lapse function a at different times. 



for the regularized case. The system can now be evolved 
for very long times and the origin remains well behaved. 
The plots show the evolution up to t = 12. 



B. System II 

Our second example is again a hyperbolic system built 
from the ADM equations and the BM slicing condi- 
tion, but now instead of using Da and Ka as funda- 
mental variables, we have used D = Da — 2D B and 
K = trK := Ka + 2Kb- In order to make the system hy- 
perbolic we remove the terms proportional to d r D B and 
d r K B from the evolution equations of K and D using the 
constraints. We then obtain a new strongly hyperbolic 
system different from System I above. 

For the regularization procedure we introduce again 



the variable A. Since in this case we have used the mo- 
mentum constraint to modify the evolution equation for 
D, we need to replace this variable with 



U :=D- 



4BX 



The final set of dynamical variables is then 

{a,A,B,D a ,U,D B ,K,K B ,X} , 

and their evolution equations are: 

d t a = -a 2 fK , 
d t A = 2aA{2K B - K) , 
d t B = -2aBK B , 
3 t D a = -8 r (afK), 
d t U = -28 r (aK) + 4aD B (K - 3K B ) 

+ 8a 

d t D B = -2d r {aK B ) 



D a K B + ^(3K B -K) 



d,K 



a 



K B + 6K% 
AXB" 



2D n 



+ K 2 



d t K t 



+ 
+ 



a 

Ar 

a 
A 
£b 
4 
2aA 
~~B~ 



U 
~2 



A 

2XB 
~A~ 



- Di 



Ar 

£l _ d r D a 
A A 

D r 



(4.18) 

(4.19) 

(4.20) 
(4.21) 
(4.22) 
(4.23) 

(4.24) 
(4.25) 

,(4.26) 



A 



D a Di 



d r D T 



U+^pj +A KK B 



d r K f 



^(K-3K B ) 



(4.27) 
(4.28) 



In this case there is no need to show any plots, as the 
numerical evolution behaves exactly in the way it did for 
System I, with the origin remaining well behaved. 



V. DISCUSSION 

Lack of regularity of geometric variables at the origin 
is often a problem for spherically symmetric evolution 
codes in numerical relativity. We have shown that the 
problem can be traced to the existence of two types of 
regularity conditions at the origin. In the first place, 
there are regularity conditions that guarantee that the 
variables are well defined at the origin. These conditions 
can be written down as a series of symmetry conditions 
at the origin for the different variables, and can be easily 
enforced in numerical simulations. However, there also 
exist regularity conditions related to the condition that 
spacetime must be locally flat at the origin. Together, 
all these regularity conditions imply that we have more 
conditions to satisfy at r — than dynamical variables, 
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which means that numerically some regularity conditions 
will inevitably be violated. We have presented a generic 
regularization algorithm that is based on the introduc- 
tion of an auxiliary variable that absorbs the problematic 
terms and on which we can impose the extra boundary 
conditions at r = 0. Our algorithm is similar in spirit, 



if not in detail, to an algorithm presented by Arbona 
and Bona for the particular case of the Bona-Masso for- 
mulation Q. We have also shown the effectiveness of 
our algorithm on a couple of different formulations of the 
evolution equations. 
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